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Abstract The orbital evolution of a dust particle under the action of a fast interstellar 
gas flow is investigated. The secular time derivatives of Keplerian orbital elements and 
the radial, transversal, and normal components of the gas flow velocity vector at the 
pericentre of the particle's orbit are derived. The secular time derivatives of the semi- 
major axis, eccentricity, and of the radial, transversal, and normal components of the 
gas flow velocity vector at the pericentre of the particle's orbit constitute a system of 
equations that determines the evolution of the particle's orbit in space with respect to 
the gas flow velocity vector. This system of differential equations can be easily solved 
analytically. From the solution of the system we found the evolution of the Keplerian 
orbital elements in the special case when the orbital elements are determined with 
respect to a plane perpendicular to the gas flow velocity vector. Transformation of the 
Keplerian orbital elements determined for this special case into orbital elements deter- 
mined with respect to an arbitrary oriented plane is presented. The orbital elements 
of the dust particle change periodically with a constant oscillation period or remain 
constant. Planar, perpendicular and stationary solutions are discussed. 

The applicability of this solution in the Solar system is also investigated. We con- 
sider icy particles with radii from 1 to 10 ^im. The presented solution is valid for these 
particles in orbits with semi-major axes from 200 to 3000 AU and eccentricities smaller 
than 0.8, approximately. The oscillation periods for these orbits range from 10^ to 2 x 
10^ years, approximately. 
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1 Introduction 



Since |Poynting| ( |1903[ ) and | Robertson] ( |1937[ ) , the dynamics of dust grains has been 
investigated in many papers and its analysis is an inseparable part of astrophysics. 
The influence of the electromagnetic radiation of a central star is usually taken into 
account in the form of the Poynting- Robertson (P-R) effect ( Poynting|1903 Robertson 
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1937 Wyatt and Whipple 1950 Burns et al. 1979 Klacka 20081. Beside the P-R 



effect, the stellar wind (corpuscular radiation of the star) can also affect the motion 
of dust particles. Covariant derivation of the acceleration caused by the stellar wind 



and its effects on the dynamics of dust grains were presented in Klacka et al. ( 2009b I 



The Solar magnetic field can affect the motion of charged dust particles in the Solar 
system ( |Parker|1958| [Kimura and Mann|1998 l. Planets can capture the dust particles 



into the mean-motion resonances (Dermott et al. 1994 Reach et al. 1995'). Because 



of the relative motion of the Sun with respect to the local interstellar medium, atoms 
of the interstellar medium approach the Sun. The direction of the approach is given 
by the actual velocity of the Sun with respect to the local interstellar medium. These 
approaching atoms form an interstellar gas flow. This flow of interstellar atoms through 



the Solar system has been already investigated in the past (e.g. Fahr||1996 Mobius et 



al.|2009 Alouani-Bibi et al.|2011 \. This interstellar gas flow can affect the dynamics of 



dust grains in outer parts of the Solar system. The assumption that dust particles in 
orbits around other stars can be also affected by interstellar gas flow is supported by 
the recent detection of debris disks around stars with asymmetric morphology caused 
by a fast motion of the disk through a cloud of interstellar matter ( Hines et al.||200T 
|Debes et al"!][2009l ) . 

The influence of interstellar gas flow on the dynamics of a spherical dust particle 
was investigated by Scherer (20001, who calculated the secular time derivatives of 
the particle's angular momentum and the Laplace-Runge-Lenz vector caused by the 
interstellar gas flow. He has come to the conclusion that the particle's semi-major axis 
can increase exponentially (Scherer 2000 p. 334). This result contradicts results of 



Pastor et al. (2011 1. Here the secular time derivatives of the Keplerian orbital elements 



of the dust particle under the action of a fast interstellar gas flow were for the first 
time calculated for arbitrary orientations of the orbit. Pastor et al. (2011 1 states that 
the secular semi-major axis of the dust particle must decrease under the action of 
the interstellar gas flow. The rate of decrease is proportional to the semi-major axis. 
Decrease of the semi-major axis only happens in secular time-scales. Therefore, the 
exact calculation of the secular time derivative of the semi-major axis was necessary. 
This result was conflrmed also by Belyaev and Rafikov (20101 who investigated the 
motion of a dust particle in the outer region of the Solar system behind the solar wind 
termination shock. They calculated the orbital evolution of the dust particle under 
the action of a constant mono-directional force, i.e., they solved a case of the classical 



Stark problem (for more information about the Stark problem see Lantoine and Russell 



|2011| and the references therein). In the Stark problem it is assumed that the orbital 
speed of the dust grain with respect to a central object can be neglected in comparison 
with the speed of the interstellar gas flow. The relative velocity terms are in |Pastor| 
et al.| (20111 taken into account to the flrst order of accuracy in the calculation of 



the secular time derivatives of Keplerian orbital elements. [Belyaev and Raflkov| ( f2010[ ) 
reproduced the result of Pastor et al. ( 2011 1 on the secular evolution of the semi-major 
axis of the particle's orbit. However, the case of a mono-directional force (i.e., the Stark 
problem) can be important for those stars where the relative speed of the neutral gas 
with respect to the central star is high. Such a situation can occur, for example, in 
merging galaxies when a star from the first galaxy moves through a molecular cloud of 
the second galaxy. 

[Belyaev and Rafikov] ( |2010| ) used the orbit-averaged Hamiltonian method. In this 
paper we use a method based on orbit-averaged time derivatives of Keplerian orbital 
elements. Using this different method we confirm results of Belyaev and Rafikov ( 2010 1. 
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We not only use a new method, we also find new results. We find an explicit form for 
the time dependence of all Keplerian orbital elements determined with respect to a ref- 
erence plane perpendicular to the gas fiow velocity vector. Mainly the time dependence 
of the longitude of the ascending node will represent a generalization of |Belyaev and| 



Rafikov ( 2010 1 results. We find a transformation for the orbital elements determined 



with respect to the reference plane perpendicular to the gas flow velocity vector into an 
arbitrarily oriented reference plane. We determine the maximal and minimal values of 
eccentricity for numerous special cases. We study in detail the behaviour of the solution 
for the planar case and for the case when the gas flow velocity vector is perpendicular 
to the line of apsides. Properties of the solution in the Solar system are discussed. 



2 Equation of motion 

In order to find the acceleration of a spherical dust particle caused by an interstellar gas 
fiow we will assume (a) that the dimensions of the dust particle are small in comparison 
with the mean free path of the interstellar gas atoms, (b) that the mass of the dust 
particle is large in comparison with the mass of the interstellar gas atoms, and (c) that 
the molecules are specularly or diffusely reflected from the surface of the dust particle. 
These assumptions lead to the following acceleration ( Baines et al.|1965 I 



^ = - C£, 7// |v - vh\ (v - vh) , (1) 

where v is the velocity of the dust grain in the stationary frame, v// is the constant 
velocity vector of atoms in the flow in the stationary frame, is the drag coefficient, 
and 7// is the collision parameter. For the collision parameter we can write 

IH = nn A , (2) 

m 

where m[{ is the mass of the neutral hydrogen atom, n/f is the concentration of the 
interstellar neutral hydrogen atoms, and A = tvR^ is the geometrical cross section 
of a spherical dust grain of radius R and mass m. Assumption (a) is a reasonable 
approximation to conditions in interstellar space. Therefore, Eq. |l| is valid for all 
dust particles with dimensions for which assumption (b) holds. 

In order to find the final equation of motion we further assume (d) that the dust 
particle orbits around a star. Therefore we take into account also the electromagnetic 
radiation of the star in the form of the P-R effect. The equation of motion then has 
the form 

dv /i 

at 

LV c / c 

-cz5 7h|v-vh| (v-v//) . (3) 



where n = GM, G is the gravitational constant, M is the mass of the star, r is the 
position vector of the particle with respect to the star, r = jrj, eji — r/r, and c is the 
speed of light in vacuum. Parameter jS is defined as the ratio of the electromagnetic 
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radiation pressure force and the gravitational force between the star and the particle 
at rest with respect to the star 



P- 



3 L 

16 TV c fi R g 



(4) 



Here, L is the star luminosity, Qp,. is the dimensionless efficiency factor for radiation 
pressure integrated over the star spectrum and calculated for the radial direction (Qpr 
= 1 for a perfectly absorbing sphere), and g is the mass density of the particle. The 
second term represents the P-R effect neglecting terms of higher orders in v/c than the 
first (|Klacka|2008[ ). Eq.|3] was used in numerical experiments by| Marzari and Thebault| 
(20Tir 



We also restrict the possible speeds of the interstellar gas flow. We will assume (e) 
that the speed of the interstellar gas flow is much greater that the mean thermal speed 
of the gas in the flow. For interstellar gas flow speeds Vf{ — |v//| comparable to the 
mean thermal speed of the gas cjj is & function of vh ( Baines et al.|[l965 |. However, 
cu has an approximately constant value for those interstellar gas flow speeds much 
greater that the mean thermal speed of the gas ( Baines et al.|[l965 Banaszkicwicz et| 
al.|1994[ [Scherer|2000[ [Klacka et al.|2009a[ ). 

We also assume (f) that the speed of the interstellar gas flow is much greater than 
the speed of the dust grain in the stationary frame associated with the central star. 

Finally, we assume (g) that the secular time derivatives of the semi-major axis 
and the eccentricity of the particle's orbit caused by the P-R effect during orbital 
evolution are low in comparison with the values of the semi-major axis and the ec- 
centricity, respectively. This assumption is reasonable for grains with large semi-major 



axes and small eccentricities (Wyatt and Whipple 19501. These particles have small 



orbital speeds. Therefore, in the P-R effect, the terms depending on velocity can be 
neglected in comparison with the Keplerian term. 

Assumptions (a)-(g) lead to the final equation of motion 



dv fi (1- P) 

^ ef( + CD IH VH Vjy 



dt 



9 e_R 



(5) 



If the radial stellar wind is also considered in Eq. (|3|, then Eq. ([sjl remain unchanged 
( Klacka et al.|[2009b| . However, assumptions (f) and (g) require radial distances that 
are probably larger than the radial distance to the stellar wind termination shock for 
the majority of Sun-like stars. 



3 Secular evolution of Keplerian orbital elements 

We want to find the influence of the constant acceleration given by the second term in 
Eq. ([5| on the secular evolution of a particle's orbit. We will assume (h) that constant 
acceleration can be considered as a perturbation acceleration to the central acceleration 
caused by the gravity of the central star. For the applicability of the assumptions (f )-(h) 
in the Solar system we refer the reader to Appendix |A] We denote the components of 
the hydrogen gas velocity vector in the stationary Cartesian frame associated with the 
central star as vh ~ (vhXj ^hYj ^hz)- Iii order to compute the secular time derivatives 
of the Keplerian orbital elements (a, the semi-major axis; e, the eccentricity; lu, the 
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argument of the pericentre; fi, the longitude of the ascending node; i, the inclination) 
we want to use the Gauss perturbation equations of celestial mechanics. To do this, we 
need to know the radial, transversal, and normal components of acceleration given by 
the second term in Eq. (|5|. The orthogonal radial, transversal, and normal unit vectors 
of the particle on a Keplerian orbit are (e.g., |Pastor|2009| | 



ei?, = 



Sat 



(cos cos(/ + oj) — sin i7 sin(/ + a;) cosi , 
sini7 cos{f + iu) + cos n sin(/ + oj) cosi , 
sin(/ + uj) sini) , 

(— cosi7 sin(/ + oj) — sin ^2 cos(/ + i^) cosi 
— sin n sin(/ + iu) + cos H cos(/ + cj) cos i , 
cos{f + u}) sini) , 
(sinJ7 sini, — cos J7 sini, cosi) , 



(6) 



(7) 
(8) 



where / is the true anomaly. For the radial, transversal and normal components of the 
perturbation acceleration, we obtain 



a-R = a vh ■ eR = a A , 
ax = a Wfj ■ ej' = a B , 
ajv = Q '^H ' = ct C 



(9) 
(10) 

(11) 



Now we can use the Gauss perturbation equations of celestial mechanics to compute 
the time derivatives of the orbital elements. The perturbation equations have the form 
(cf., e.g., [Murray and Dermott|1999| [Danby||1988 1 

[a^ e sin/ + ay (1 + ecos/)] , 



da 


2 a 


dt ^ 


l-e2 


de 




di ^ 


V 


doj 


1 


~di ~ 


e V 



a^ sin / + ay I cos / 



e + cos / 



dn 

~dt 
di 
dt 



a cos / — ax 
sin(/ + cj) 



1 + e cos / 
h e cos / 



1 + e cos / 



sin / 



sm I 
sin(/ + oj) 



r 



ajV cos(/ + uj) , 



where p = a(l — e ). The time average of any quantity 
can be computed using 

{9) 



(12) 

during one orbital period T 



df 



2 TT a3/2 
1 



df 



2tv a'^ VT 



e2 Jo 



g r df 



(13) 
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Fig. 1 A schematic representation of the values S, I and C for a given orbit. 



where the second {^/Jtpp = r df /dt) and the third (47r a — ^pT ) of Kepler's laws 
were used. From Eqs. (|9|-( 13 1 we finally obtain for the secular time derivatives of the 
Keplerian orbital elements 



(14) 
(15) 
(16) 
(17) 
(18) 




where the quantities 
5" = 



C 



(cos Q cos Lo — sin Q sin lo cos i) vjjx 
+ (sin n cos uj + cos H sin to cos i) v^y 
+ sintj sini V}{z i 

(— cos J7 sinoj — sinJ? coscj cosi) vjjx 
+ (— sini7 sincj + cosf? cos a; cosi) v^jy 
+ cos UJ sin i vjjz > 

sin il sin i vux ~ cos Q sin i Vfjy + cos i vu z 



(19) 



are values of yl = ■ e/j, B — V[{ ■ ex and C — wu • e^y at the pericentre of the 
particle orbit (/ = 0), respectively. The value of C is a constant on a given oscular 
orbit. The values of S, I and C are depicted in Fig. [l] For the special case v// = 
(—V[j, 0, 0) the system of Eqs. ( 14 )-{ 18 1 and Eqs. ( |19[ ) reduces to the system of Eqs. 
(3)-(7) in Mignard and Henon| ( 1984) with the exception of one term in the secular 
time derivative of the longitude of ascending node. iMignard and Henon ( 1984 1 studied 



the motion of a particle around a planet subject to radiation pressure from a central 
star. Mathematically, they considered the Stark problem in a rotating reference frame. 
The extra term in the secular time derivative of the longitude of ascending node is 
caused by the Coriolis force. The Stark problem is recovered if the rotation rate is 
set to zero. Because of the extra term caused by the Coriolis force their analytical 
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solution is different and does not correspond to the solution which will be given in this 
paper. Their solution of the Stark problem is not easy to see from their notation in the 
rotating reference frame. Moreover, they considered only the case for which v// lies in 
the plane with i — and is antiparallel with x axis. Our solution will be more general. 



4 Secular evolution of orbit with respect to gas flow velocity vector 



The parameters S, I and C determine the position of the orbit with respect to the 
gas flow velocity vector. Therefore, their time derivatives are useful for a description of 



the evolution of the orbit's position in space. Putting Eqs. ( 16 1-( 18 1 into the formulas 
for the averaged time derivatives of the quantities S, I and C given by Eqs. ( |19[ ), we 
obtain 



dt 

dl 
dt 

dC 
dt 




e 



(20) 
(21) 
(22) 



Eqs. (|20|-(|22| are not independent, because S (dS/dt) + I (dl/dt) + C (dC/dt) = 
always hold. Eqs. (20l-(22l together with Eqs. (14l-(15l represent the system of 



equations that determines the evolution of the particle's orbit in space with respect 
to the gas flow velocity vector. All orbits that are created with rotations of one orbit 
around the line aligned with the gas flow velocity vector and going through the centre 
of gravity will undergo the same evolution determined by this system of equations. 

Now we find the solution of the system of equations given by Eqs. ([T4|-([l5| and 
Eqs. (I20|-(l22|. If we combine Eq. (fT5| with Eq. pO|, then we obtain 



dS 
dt 



de\ S 
dt , 



This equation leads, for 5 7^ 0, to the differential equation 

de 



dS 



with the solution 



D 



(23) 



(24) 



(25) 



Here, _D is a constant which can be determined from the initial conditions. Thus, if 
the major axis of the orbit is aligned with the direction of the gas flow velocity vector, 
then the eccentricity is minimal. If we combine Eq. (|15[) with Eq. (|22[), then we obtain 



dC 
~dt 



e C 



This equation leads, for C 7^ 0, to the differential equation 

dC e de 



(26) 



(27) 
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with the solution 



\C\ 



(28) 



where F is an integration constant. Thus, if the magnitude of the normal component 
of v/f measured in the perihelion is maximal, then the eccentricity is maximal. For / 
we obtain from S" 



+ 7^ + 



'H 



„2 



F2 



(29) 



If S = 0, then from Eq. (201 we obtain {dS/dt} = 0. Since Eq. (281 holds for the case 



5 = and C / 0, also Eq. (29 1 with D = holds. Therefore we put D = in Eq. 



(25} for the case S = and C / 0. If C = 0, then from Eq. (|22| we obtain (dC/dt) 
= 0. Since Eq. (|25| holds for the case C = and 5* / 0, also Eq. ([29| with F = 
holds. Therefore we put F 

C = 0, then we have |/| — vh- Therefore we put D 
Eq. (28 1, respectively. We come to the conclusion that Eqs. (25 1, (28 1 and (291 always 

Eq ~ 



in Eq. ([28| for the case C = and S* 7^ 0. If S' = and 

and _F = in Eq. pSl and 



together with 



hold during the orbital motion of the particle. Eq. ( |25[ ) and Eq. (28 
the properties of the system of differential equations given by Eqs. ( 14 
( 20 1-( 22 1 imply that S and C can not change sign during the orbital evolution of the 



( 15 1 and Eqs. 



particle. If _D 7^ 0, then both S and e must be non-zero during the orbital evolution of 
the particle. Similarly, if F 7^ 0, then C 7^ and e 7^ 1 during the orbital evolution of 
the particle. For the special cases D = or F = is necessary to use the properties of 



the whole system the differential equations given by Eqs. ( 14 1-( 15 1 and Eqs. ( 20 1-( 22 \ 
because, as we will see later, in these cases the eccentricity can be close to or 1 
respectively. Therefore we can write 



5 - 



U 
e 



(30) 



C 



V 



and 



|/| = 



^2 

g2 i__ g2 



(31) 



(32) 



where U and V are some constants. We come to the conclusion that also Eqs. (|30| 
always hold during the orbital motion of the particle. 

The S, and C components of the hydrogen gas velocity vector at the pericentre 
of the particle's orbit depend only on the particle's eccentricity. To find the evolution 
of the eccentricity, we insert Eq. (1321) into Eq. (1151). We obtain 



de 
dt 



= ± 



3 a V}j 
2 e 



;2(1- 



V2 



(33) 



The plus sign is for positive values of / and the minus sign for negative values of /. 
If I is negative, then the eccentricity decreases. If I is positive, then the eccentricity 
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increases. In order to integrate this equation we rewrite the expression in the second 
square root in the following form 



— 1 - e 



{el - e')(e^ - el) 



where 



2 

ei 



1 + 



+ 



1 + 



(72 y2 



t/2 



2 

62 = 



1 + 



V2 



(72 



(34) 



(35) 



(36) 



It is possible to show (using Eqs. 30 31 and 32 1 that the expres sion in the square root 
is always positive or zero, el G [0, 1] and e| G [0, 1]. From Eq. |33l we obtain for / ^ 




/ 



e de 

^(e2_e2)(e2-e|y 



= ± 



3 a V}j 



"2 
a 



(37) 



where is an integration constant. Hence 



2,2 2 2 
ei +62 ei - 62 



cos I ± ZaVfj 



(38) 



From this equation we can see that the eccentricity changes between its minimal value 
62 and its maximal value 6i for one solution determined by the choice of the sign (+ 
or — ). We denote the time close to the maximum eccentricity ei as tjv/- From Eq. (38 1 
we obtain for values of ip at the time 



2((5+ = 3 a Vf{ 



i-M — vr — 2 fci TT , 



(39) 



2(/?_ = — 3 a Vff 



— vr — 2 ^2 TT 



(40) 



where ki and k2 arc two integers. In Eq. (38 1 both these values lead to the same 
solution 

9.9 9 9 / r 

a 



9 e? + en 6? — 6o 
e = Z + „ cos 3avH 



{t - tM) 



(41) 



This solution is in accordance with the result of Belyaev and Rafikov (20101 up to the 
definition of an additive constant for the time. The eccentricity changes periodically 
with the oscillation period 

(42) 

3 a vh \ a 

For a connection between this result (obtained from Eq. |5| and the results obtained 
from the equation of motion with the relative velocity v included in the force caused 
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R = 


= 1 fim 


R = 


= 2 fim 


a 


Te 




a 


T 


[AU] 


r 10^ 1 

[years J 




[AU] 


r lo'^ 1 

[years J 


200 


2.06 




200 


5.35 


300 


1.68 




400 


3.78 


400 


1.46 




600 


3.09 


500 


1.30 




800 


2.67 


600 


1.19 




1000 


2.39 


700 


1.10 




1200 


2.18 



R = 


= 5 fj,m 


R = 


10 /^m 


a 


Te 




a 


T 


[AU] 


[years J 




[AU] 


r lo'^ 1 

[years J 


200 


14.90 




500 


19.45 


600 


8.60 




1000 


13.75 


1000 


6.66 




1500 


11.23 


1400 


5.63 




2000 


9.73 


1800 


4.97 




2500 


8.70 


2200 


4.49 




3000 


7.94 



Table 1 Oscillation periods Te determined for orbits with various semi-major axes in the 
Solar system. Dust particles with radius i? £ {1 fim, 2 /im, 5 fim, 10 /^m}, mass density q = 
1 g cm~^ and Q' = 1 are used. 



by the interstellar gas flow, we refer the reader to Appendix |B] The values of Te for 
dust particles in the Solar syste m are in Ta ble [l] For the Solar sy stem we used 



^26 



(e.g. Lallement||1996 Landgraf ct al 



20021, cj) = 2.6 ( Banaszkiewicz et al. 



1999 1 (see Appendix 



1994 



Scherer 



2000 



Pastor et al.|2011|, nn = 0.2 cm~^ (iBelyaev and Rafikov|2010|) and vh = 26 km s 
_ i ^ 



To find the time evolution of S, \I\ and C, we can put Eq. (41 1 into Eqs. (30l-(32l. 
S, \I\ and C also change periodically with period Te. Now we find the evolution of / 
from Eq. (|32|. If we put Eq. (|4l| into Eq. (l34|, then we get 



{el - e^)(e^ - el) 



4 



X sin^ \ 3avH J— (i — *m) 
V 



(43) 



Hence, 



VH 



2 2 
61-62 



sin I 3av[{ 



(t - tM) 



(44) 



If I is negative, then the secular eccentricity must decrease (see Eq. 15 I. Therefore if 



we compare the evolutions given by Eq. (41 1 and Eq. (441, we come to the conclusion 
that 



I = 



62 



sin I Zavjj 



(45) 



Eq. (41 1 and Eq. (45 1 hold also for the special cases 17 = or F = 0. 



4.1 Planar case 



C = for the special case when the velocity of the hydrogen gas vjy lies in the orbital 
plane of the particle. Putting C = (V^ = 0) in Eqs. (35l-(36l, one gets 



2 
ei 



(46) 



2 



(47) 
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Therefore the minimum eccentricity is \U\/vu and the maximum eccentricity is 1. 

(48) 



Using Eqs. (46l-(47l we obtain from Eqs. (|30|, (41 1, and (45l, for the planar case 



vh 



— sin I Savff 



2 1 

e =1 



1 - 




(t-tM) 



(49) 
(50) 



We can mention that S'^ + 7^ = v'j^ always holds for the planar case. Eq. (501 is in 



accordance with the result from Belyaev and Rafikov (20101 for the planar case up to 



the definition of an additive constant for the time. 



For the planar case we obtain from Eq. (21 1 



> 



(51) 



M/3 e 

Therefore, in the planar case I increases with time. This is in accordance with Eq. 



( 49 1 . To show this we rewrite Eq. (|49|l as 
/ 



e 



(72 



cos I 
4 V 2 





1 a 


2 \ 






/ a 


2 \ 





(52) 



If the eccentricity e is close to ei = 1 (see Eq. 



to 



50 1 , then I changes from 



vy — U"^. Because of this property, / can always increase. Thus, in the planar 



case the orbit rotates into position with a maximal value of I. As the eccentricity 
increases, the dust particles gets closer to the central star, because the semi-major axis 
is constant. We must note that this theory is less applicable for larger eccentricities 
(see Appendix [A|). 



5 Secular evolution of Keplerian orbital elements determined with respect 
to a plane perpendicular to gas flow velocity vector 

The choice of coordinate system was up to now arbitrary. We will denote with two 
primes those coordinate systems in which the gas flow velocity has a vector direction 
aligned with the direction of the z"-axis. In such a coordinate system we have vjj = 



{0,Q,vh). From Eqs. (191 we obtain 



r, . a ■ .11 

a = smtj smz vj^ , 

J II ■ -11 

1 = cos cj sm I Vfj , 

C = cosi" Vh . (53) 
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From these equations we immediately obtain 



tana; 



S 
" / ' 

C_ 

VH 



(54) 



(55) 



The last unknown orbital element in this coordinate system is O" . If we use Eqs. (531 



in Eq. ( 171, then we can obtain 



dn" 

dt 



3 a 



S C 



(56) 



2 fff V sin^ i" 1 — ' 
From the equation above we can deduce that dQ" /dt is always positive, negative, or 



zero because S and C can not change sign (see Eq. 30 and Eq. 31 1 . To find the evolution 
of fl" , we divide Eq. (56 1 by Eq. ( 15 1. If we use in the result of division Eqs. ( 30 1, ( 31 1, 



(34 1 and (551, then for / 7^ we finally obtain 
dQ" U V e 



vl ^(ef-e2)(e2-ei) 1 - 



(57) 



The minus sign is for positive values of / and the plus sign is for negative values of /. 
Integration of this equation yields 



± il = TjyyT ^.rctan 



y2 



^2 ef 



y2 



+ 1/' 



(58) 



Now the plus sign is for positive value s of / and t he minus sign is for negative values 
of /. It is possible to show (using Eqs. 



30 



31 



and 



32| that 1 - V'^/vjj 



Ci > and 

1 — V^/vfj — 62 — 0. If we insert Eq. (411 into Eq. (581, then, after some algebraic 
manipulations, we finally obtain 



tan 



\UV\ 

UV 



± n" 



A/2 



V2 



tan(^ — 



{t-tM) 



(59) 



This equation leads to two equations (compare Eqs. 45 and 59 1. One for J > 
tan (f?" — : 



U V 
\UV\ 



y2 



( 3avH 
X tan I — ^ 



(60) 
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and one for J < 



tan ( — ^7 — tp-j = 



U V 

\UV\ 



V 



( SavH 
X tan I — ^ 



{t - tM) 



(61) 



These solutions meet at time when I changes its sign (see Eq. 45 1. Close to the 
time tM, the right-hand side of both equations is close to zero. For values of and 
ip- at the time £m we obtain 



M 



— ^4 TT 



(62) 

(63) 



where ^3 and are two integers and O'^ is the value of Q" close to the time tj^j- In 



Eqs. (60 1 and (61 1 both these values lead to the same solution 



tan (^n" — fS'l.i^ 



U V 



\uv\ 



y2 



X tan 



y2 



(64) 



This equation is a generalization of the result of Belyaev and Rafikov (20101. From the 
discussion of Eq. ( 56 1 we know that fi" is a monotonic or constant (for 5* = or C = 



and e 7^ or e 7^ 1) function of time. If we consider values of O" only in the interval 
[0, 2n), then n" is a periodic or constant function of time. 



5.1 Stationary solution 



Eqs. (20l-(22l and (30l-(32l enable finding a stationary solution determined by two 

(65) 



equations 
and 



= 

e 







(66) 



If we insert Eqs. ( |30| and (|31| into Eq. (|66|, we obtain a condition for the eccentricity 
of the stationary solution 

2 \U\ 



\u\ + \v\ 



(67) 



If Eqs. (651 and (67 1 are fulfilled, then e, 5*, / and C remain constant. 

For orbital elements determined with respect to the plane perpendicular to the gas 
fiow velocity vector, we obtain from condition / = in the second of Eqs. ( 53 1 that 
cos a;" = (since sini" 7^ 0). Hence 



UJ = - + fcsTT 



(68) 
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Fig. 2 Transformation from a frame in which the 2:"-axis is ahgned with the interstellar gas 
velocity vector into an arbitrary reference frame. 



where fc5 is an integer. Eq. ( 18 1 yields that the inclination is in this case a constant. 
Finally from Eq. ( 17 1 we obtain that fl" depend linearly on time, fl" is the only orbital 
element of the stationary solution that varies with time. Therefore, the orbital plane 
rotates around a line aligned with the gas flow velocity vector and going through the 
centre of gravity. This stationary solution is in accordance with the result of IBelyaev] 



and Rafikov (20101. 



6 Evolution of orbital elements determined with respect to an arbitrary 
reference plane 

It is useful to have the evolution of orbital elements in a reference frame oriented 
arbitrarily with respect to the hydrogen gas velocity vector. Therefore in this section 
we transform the orbital elements derived in the previous section into such a frame. 

The frame denoted with two primes has the z"-axis aligned with the velocity vec- 
tor of the hydrogen gas. x" is the axis from which fl" is measured. For the sake of 
simplicity we will assume that the x"-axis lies in the a;i/-plane of the frame into which 
we want to transform the orbital elements. By making this assumption we do not lose 
any generality of the solved problem since the direction of the a;"-axis is arbitrary. 
The transformation of the particle's coordinates form the two primed frame into the 
unprimed frame can be made by the composition of two rotations (see Fig. |2|. The 
first rotation is around the a::"-axis by the angle between z' and z" . 

y = y cos (pr + z sm (pr , 

z = — y" sin <f)r + z" cos (f>r , (69) 

where 

sm(j)r — — , 

VH 

COS 0r = — . (70) 

VH 
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The second rotation is a rotation around the z'-axis by the angle between x and x' . 

X = x' cos tpr + y sin tl)r , 
y = — x' sin t/ir + y' cos 4>r , 

z = z' , (71) 



where 

sin tpr 



i2 I T,2 



cosVr = "-^^ . (72) 



"HX ^ "HY 

Hence the transformation from the two primed frame into the unprimed frame is 

X — x" cos tpr + {y" COS ipr + z" sin cpr) sin ipr , 

y — — x" sin t/ir + {y" cos <^r + z" sin 0^) cos 'ipr , 

z — — y" sin <jf>r + cos 0r , (73) 

In order to find fi and i we can transform the coordinates of the unit vector = 
(sin J?" sini", — cosi7" sini", cosi"). We obtain 

ejvx — sm U sm i 



HX ' "HY 



+ I — cos U sm I h cos J — 

I VH VH 

X """^ , 

\J ^HX + ^HY 

■ o" • ■" ^HX 

ej\[Y ~ —smU smi 



+ I — cos U sm I h cos I — 

I VH VH 

VHY 



2 I ,,2 
HX ^ "HY 



ejvz ~ cos J / sm i h cos j . (74) 

VH VH 



In this equation, Eq. (64 1 and Eq. 1 55 1 have to be used to determine the values of J?" 
and i" . Since ejy = (sinf? sini, — cos J7 sini, cosi), we can calculate O and i from 

tan ^2 = - , (75) 

eNY 

cosi = ejvz • (76) 
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Now we can find oj from Eqs. ( 19 1. One can easily verify that 



Therefore 



J — S cos u — I sin ij = cos il v^x + sin J7 v^y i (77) 

H — S sin (jj + I cos (jj — — sin f2 cos i 

+ cos O cos j Vhy 

+ sin i . (78) 



HS-JI , JS + HI 



6.1 Perpendicular case 



If the velocity vector of a neutral gas is perpendicular to the line of apsides, then 5 
= ([/ = 0) and the value of S does not change with time. For this special case we 
obtain from Eqs. ( 35 1-( 36 1 



2 1 

ei = 1 



2 n 
62 = 



(80) 



(81) 



Therefore the minimum eccentricity is and the maximum eccentricity is y^l — V'^/v'j^. 
Using Eqs. ([80|)-(|8T| we obtain from Eqs. ([3T|, (|4T| and ((45|) for the perpendicular case 



C - 



V 



71 ' 



(82) 



e2 2 

2 



— sin ( Savff 



(t - tM) 



V 



H 



I Savff r~a 



{t - tM) 



(83) 
(84) 



We mention that + l"^ = vfj always holds for the perpendicular case. 



For the inclination in the two primed frame, we have from Eq. (31 1 and Eq. (55 1 



cos J — 



VH 



V 

VT- 



(85) 



Hence for i" we obtain that i" £ (0, a.TCCos{V / v h)) for positive values of V and i" £ 
(arccos(V/f/f ), vr) for negative values of V . 

Since S = and sini" 7^ 0, we obtain from the first of Eqs. (53 1 that sincj" = 0. 
Therefore 

Lu" — fcgTT , (86) 

where fcg is an integer. 

Finally Eq. (561 for sini" / yields that O" is a constant. 
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Now we find relations between the orbital elements in the unprimed frame for the 

(87) 



special case 5* = 0. If we combine Eq. ( 16 1 for 5 = with Eq. ( 18 1, we get 

cos UJ 



cos I 



di 



which has the solution 



smi smco] = 



(88) 



where d is an integration constant. Since i £ (0, vr), sini > 0. If sinoj 7^ initially, then 
d =^ and therefore sin a; can not change sign. Hence 



Sim smoj : 



I, 



(89) 



where I is a constant. From Eq. (89 1 we can see that the orbital plane is only tilted 
back and forth around the line of apsides for the special case S — 0. The instantaneous 
value of i can be determined from Eq. (76 1. In order to find H we combine Eq. (17 1 
with Eq. (|18[). We get differential equation 



coscj aSJ = at 



(90) 



Equation (89 1 allows us to convert this expression into an equation for with 
solution 

,^ ^ , sincj cosi ,„,^ 

tan(^? - Oq)^ , (91) 



where Hq is the value of i? for i — ■k/2. 



7 Conclusion 

We have investigated the long-term orbital evolution of a spherical dust particle per- 
turbed by a small constant mono-directional force caused by a fast interstellar gas flow. 
The secular time derivatives of the particle's Keplerian orbital elements were derived. 
We transformed the system of differential equations for the Keplerian orbital elements 
into a system of differential equations for the semi-major axis, eccentricity, and the 
radial, transversal, and normal components of the interstellar gas fiow velocity vector 
determined at the pericentre of the particle's orbit. In these new variables, the system 
of differential equations can be easily solved. The solution of the system was used in 
order to obtain the evolution of the Keplerian orbital elements in a reference frame in 
which the orbital elements are determined with respect to the plane perpendicular to 
the interstellar gas velocity vector. We found an explicit form for the time dependence 
of all Keplerian orbital elements in this reference frame. We generalized the expres- 
sion for the time dependence of the longitude of the ascending node found by |Belyaev| 



and Rafikov (20101. We transformed newly found orbital elements into an arbitrary 
reference frame. This transformation gave us explicit time dependences of all Keple- 
rian orbital elements in the arbitrary reference frame. The orbital elements of the dust 
particle in an arbitrary reference frame change periodically with a constant oscillation 
period or else remain constant. We determined the properties of the solution for the 
planar case and for the case in which the gas flow velocity vector is perpendicular to the 
line of apsides. We found the maximal and minimal values of the eccentricity in these 
cases. In the planar case, the particle's orbit approaches the direction with maximal 
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1 2 fim 




«[■] i>[-l 



Fig. 3 Allowed values of orbital semi-major axes and eccentricities in the ae-phase plane 
dete r mine d usi ng t he relations between the semi-major axis and eccentricity obtained from 
( |93| , \95\ and l |96| . Dust particles have radius ij £ {1 nm, 2 /xm, 5 ti m, 1 /xm} and mass 
den sity (? = 1 g cm"'' and Qp^ = 1. Dashed line indicates |93|, solid line, |95l and dotted line, 
||96|. The allowed region of the semi-major axes and ecccritricities is depicted in grey. 



value of I. In the perpendicular case, the orbital plane is tilted back and forth around 
the line of apsides. We also confirmed the stationary solution found by |Belyaev and| 



Rafikov (20101. For the stationary solution, the orbital plane rotates around the line 
aligned with the gas fiow velocity vector and going through the centre of gravity. 

This solution can be applied also for the dust particles in the Solar system. If we 
consider icy particles with radii from 1 to 10 ^im, then the solution is valid for orbits 
with semi-major axis from 200 to 3000 AU, approximately. More exact values of these 
limits depend on the radius of the particle. A maximal orbit eccentricity for which the 
solution is valid is smaller than approximately 0.8 and depends on the semi-major axis 
of the orbit and the particle's radius (see Fig. [sf. The period of change of the orbital 
elements for these orbits ranges from 10^ to 2 x 10^ years, approximately. 



A Applicability of assumptions (f)— (h) in Solar system 

Assumption (h) can be mathematically written as 

fip CD riH niH A vjj 



(92) 
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This inequality must always hold during orbital evolution. The right-hand side of the inequality 
is a constant and the left-hand side is minimal at maximal r. Since we consider elliptical orbits, 
r is maximal at the apocentre of the orbi t. T herefore we will check this inequality at the 
apocentre. In the apocentre we can rewrite | |92[ l as 

3 CO nu ran vjj 
a2 (1 -h e)2 4 H £. ' ^ ' 

where a is the semi-major axis and e is the eccentricity of the dust particle's elliptical orbit. 
For a given particle this inequality determines the allowed a and e in the ae-phase plane. 
Assumption (f) can be mathematically written as 

Wh\ » M . (94) 

Because the velocity of the particle is maximal at the pericentre of the orbit, we will check 
this inequality at the pericentre. In the pericentre we can rewrite l|94[l as 



IvhI » . (95) 

V a 1 — e 

Similarly to | |93| , this inequality also determines the allowed a and e in the ae-phase plane. 

Assumption (g) can be approximately mathematically written as ^Wyatt and Whipple] 
|1950| l 

/'^°\ 1 1^1^ 2 + 3e2 
- U/PH a = ^ (l-e2)3/2 « ' (96) 

and 



\dt/pj^ a ca'^ (l-e2)3/2 

/ de \ 1 _ 5 /3 1 
'Ut/pp e - (1 - c2)i/2 « ■ 

Here, bg and b e de termine the maximal values for the relative secular time derivatives. Inequal- 
ities l96l and d97l for a given particle, ba and fee determine the allowed a and e in the ae-p hase 
plane. It is possible to show that for bg > 5fea/4, the semi-major axis determined from | |96[ l 
is always greater than the semi-major axis determined from | |97| for e £ [0, 1). Therefore we 
will check only inequality l |96[ l. Using inequality | |96| l we take into account also the eccentricity 
of the particle's orbits. This is a more general method than the method used b y Belyac v and] 
[Rafi kov (2010) who considered only circular orbits (see Eq. 26 in |Belyaev and Rafikoy 201011 . 

T o obtain the conditi ons in the Solar system w e used = 3.842 X 10^° W |Sa.hcall| 
|2002 |, en = 2.6 ( Banaszkiewicz et al.|l994| |Scherer|20 00 Pas tor ct al. 2oTT|. uh = 0.2 cm" - 
( [Belyaev and RafikOT(2010l l and vh = 26 km s~^ (e.g. |LaIlement|1996| [Laiidgrat et al.|l999| 
In Fig. |3j four panels are depicted. Each of these corresponds to a different particle radius. 
We used particles with radius H S {1 ^m, 2 /.tm, 5 /^m, 10 /.tm} and mass density ^ = 1 g 
cm~3 and Qi^ = 1. De pend enc es between the semi-major axis and eccentricity are obtained 
from ( |93[ l, l |95| and | |96| l. In | |93| we assumed that the left-hand side must be at least 10 times 
greater than the right-ha nd s ide. For the equation obtained using this inequality we used the 
dashed line in Fig. ^ In l |95| we also assumed that left-hand side must be a t le ast 10 times 
greater than the right-hand side. In Fig.jsjthis is depicted using a solid line. In | |96[ l we assumed 
that ba = 10~® year~^ and the right-hand side must also be at least 10 times greater than 
the left-hand side. In Fig. [s] | |96[ ) is depicted useing a dotted line. The grey region in Fig. [s] 
represents values of semi-major axes and eccentricities in the ae-phase plane for which all three 
inequalities hold. In Fig. [s] we can see that assumptions (f )— (h) can be applied to the Solar 
system. If also assumptions (a)-(e) hold, then Eq. can be used. Assumption (h) enables 
the use of perturbation theory. 



B Connection with previous work 



[Pastor et al.| | |2011 1 take into account also the relative velocity v of the dust particle with 



respect to the star in the force caused by the neutral interstellar gas. The particle's equation 
of motion has the following form 

-t: = - ^ '^D IH |v - vh\ (v - vh) • (98) 
dt r'^ 
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After expansion of the right-hand side of this equation into a series the linear term in the 
relative velocity is included in the calculation of the secular time derivatives of the Keplerian 
orbital elements. The result for the secular time derivative of the semi-major axis is 



da \ 
dt / 



1 + 



(99) 



where 



(100) 



From Eq. (|99j we can see that the secular time derivative of the semi-major axis is proportional 
to the value of the semi-major axis (the value of a is independent of the semi-major 

axis). The term within the square brackets (SB) in Eq. (|99| can be written as 



SB : 



[P - 5^) 



1-^ 



(/Vl-ea + S^) > 



(101) 



The secular semi-major axis is a decreasing function of time. 

Now, we find the maximal possible decrease of the secular semi-major axis. Because terms 
the multiplied by S^ and are both positive, we obtain the maximal value of SB for the 
orbit orientation characterized by C = 0. If C = 0, then S^ + I'^ = v^. Using this, the value 
of SB can be written as 



SB : 



1 - Vi - I 



(102) 



characterized by S 



v^. Hence, the maximal value of SB is 



(103) 



In order to find the behaviour of the function g{e) we can write 



de 



(104) 



de de V 

= -2e+-^>0. (105) 

Because dgi{e)/de > 0, gi{e) is an increasing function of the eccentricity. The value of gi{0) 
is 0. Therefore, gi{e) is positive for e £ (0, 1]. If gi(e) is positive, then dg(e)/de > 0. Because 
dg{e)/de > 0, the function g{e) is an increasing function of the eccentricity for e £ (0, 1]. The 
function g{e) has its maximal value for e = 1. Therefore the maximal value of SB is 

5Bmax = vjj . (106) 

Hence, the maximal possible decrease of the secular semi-major axis is 

(^) = - 4a Co -JH vjj a . (107) 

\ dt / jnax 
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R 


Ta 


[/xm] 


1 years | 


1 


46.7 


2 


93.4 


5 


233.4 


10 


466.8 



Table 2 Time intervals Ta after which the decrease of the semi-major axis is smaller than 
10% {ca = 0.9) for particles with ij £ {1 /jm, 2 /^m, 5 /xm, 10 /im} and mass density p = 1 g 
cm~^ in the Solar system. 



We define a time interval Ta during which we suppose that the semi-major axis is approx- 
imately constant. For a given value of the semi-major axis we have 

Ta 

^ '-^^ = constant = Ca ^ I . (108) 

a 

Becaus e we use the maximal possible decrease of the secular semi-major axis, the left-hand side 
of Eq. | |108| will be even closer to unity for the real secular time derivative of the semi-major 
axis. However 

Ta 3Taa v„ 



2 n Y M/3 



(109) 



(see Eq. |42| is proportional to a^/^. Therefore, for larger semi-major axes we obtain more 
eccentricity periods in the same time interval Ta- Therefore the theory with the equation of 
motion considered in the form of Eq. ([sjl is better applicable for larger semi-major axes. We 
can calculate the values of Ta from Eq. | |108| . If we use Ca = 0.9, then we obtain for dust 
particles with R £ {1 fim, 2 fim, 5 ^m, 10 fim} and mass density (? = 1 g cm^'^ in the Solar 
system (see Appendix the values of Ta shown in Table |2] 
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